# SLE human patients
library(librarian)
shelf(Seurat,
      tidyseurat,
      tidyverse,
      harmony,
      SingleR,
      ggrepel,
      ggpubr)
source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

# examine mva pathway ----------
kegg_mva <-
  c('Hmgcr','Hmgcs1','Hmgcs2','Fdps','Mvd','Idi1','Idi2', 'Pmvk','Acat1','Acat2','Mvk') |>
  str_to_upper()

# examine interested cytokines ------
key_cytokine <- c('Il6','Ccl20','Tslp','Flt3lg','Csf2','Tnf') |>
  str_to_upper()

# read mtx data ----------
safe10x <- safely(Read10X, quiet = F)

# turn warning to error
options(warn = 2)
## zheng2022 --------
zh22_dir <- list.dirs('mission/FPP/xiangya/zheng2022/',
                      full.names = T, recursive = F)

zh22_dir <- zh22_dir |>
  basename() |>
  set_names(zh22_dir, nm = _)

batch1_zh22 <- zh22_dir[1:8] |>
  Read10X(strip.suffix = T)

batch2_zh22 <- zh22_dir[9:17] |>
  Read10X(strip.suffix = T)

zh22b1_name <- zh22_dir[1:8] |> names()

batch1_zh22[1:5,1:5]

batch2_zh22[1:5,1:5]

batch1_zh22 |> glimpse()
batch2_zh22 |> glimpse()

intsc_genes <- rownames(batch1_zh22) |>
  intersect(rownames(batch2_zh22))

## new xiangya data -----
nxy_dir <- list.dirs('mission/FPP/xiangya/',
                      full.names = T, recursive = F) |>
  str_subset('SLE')

nxy_dir <- nxy_dir |>
  basename() |>
  str_c('New') |>
  set_names(nxy_dir, nm = _)

batch_nxy <- nxy_dir |>
  Read10X(strip.suffix = T)

batch_nxy |> glimpse()

intsc_genes <- rownames(batch_nxy) |>
  intersect(intsc_genes)

## kim2023 --------
kim_dir <- list.dirs('mission/FPP/xiangya_sle_scRNA/data/kim2023/',
                     full.names = T, recursive = F)

kim_dir <- kim_dir |>
  basename() |>
  str_c('kim') |>
  set_names(kim_dir, nm = _)

batch_kim <- kim_dir |>
  Read10X(strip.suffix = T)

batch_kim |> glimpse()

intsc_genes <- rownames(batch_kim) |>
  intersect(intsc_genes)

# turn to normal warn level
options(warn = 1)

# create sobj -------
sobj <- list(batch_kim, batch_nxy, batch1_zh22, batch2_zh22) |>
  map(\(x)x[intsc_genes,]) |>
  reduce(RowMergeSparseMatrices) |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj$mito.ratio <- sobj |>
  PercentageFeatureSet('^MT-')

sobj |> VlnPlot('mito.ratio', pt.size = 0)

sobj <- sobj |>
  filter(mito.ratio < 10)

sobj <- sobj |>
  mutate(sle = str_extract(orig.ident, 'HC|SLE'),
         dataset = case_when(
           str_detect(orig.ident, 'kim') ~ 'kim2023',
           str_detect(orig.ident, 'New') ~ 'xiangya_new',
           orig.ident %in% zh22b1_name ~ 'zheng2022_batch1',
           .default = 'zheng2022_batch2'
         ),
         individual = str_remove(orig.ident, 'D$|E$'))

sobj <- sobj |>
  quick_process_seurat(batch = c('orig.ident', 'dataset', 'individual'))

sobj |> DimPlot(group.by = 'dataset')

sobj |> DimPlot(group.by = 'sle')

## auto annotate cell type ---------
hpca <- celldex::HumanPrimaryCellAtlasData()

sobj <- sobj |>
  mark_cell_type_singler(ref = hpca, new_label = 'hpca_main')

sobj |> DimPlot(group.by = 'hpca_main', cols = DiscretePalette(36),
                label = T, label.box = T)

sobj |>
  ggplot(aes(tissue, fill = hpca_main)) +
  geom_bar(position = 'fill') +
  scale_fill_manual(values = DiscretePalette(36))

sobj |>
  ggplot(aes(sle, fill = hpca_main)) +
  geom_bar(position = 'fill') +
  scale_fill_manual(values = DiscretePalette(36))

sobj <- sobj |>
  mutate(tissue = case_when(
    str_ends(orig.ident, 'D') ~ 'dermis',
    str_ends(orig.ident, 'E') ~ 'epidermis',
    .default = 'mixed'))

sobj |>
  ggplot(aes(tissue, fill = hpca_main)) +
  geom_bar(position = 'fill') +
  scale_fill_manual(values = DiscretePalette(36)) +
  theme_pubr()

# keratinocyte & melanocyte marker genes
sobj |> DotPlot(c('KRT1','KRT14','PMEL','MLANA'),cluster.idents = T) +
  coord_flip()

sobj |>
  dplyr::count(orig.ident, individual, dataset) |>
  write_tsv('human_trpv3_sle.sample.tsv')

## manual edit cell type -----
sobj <- sobj |>
  mutate(manual_main = case_when(seurat_clusters == 22 ~ 'Melanocytes',
                                 hpca_main == 'Tissue_stem_cells' ~ 'Fibroblasts',
                                 .default = hpca_main))

sobj |> DimPlot(group.by = 'manual_main', cols = DiscretePalette(36),
                label = T, label.box = T)

## rds checkpoint -------
sobj |> write_rds('mission/FPP/human_SLE_27.rds')

sobj <- read_rds('mission/FPP/human_SLE_27.rds')

kimmeta <- sobj@meta.data |> filter(dataset == 'kim2023')

kim_dc <- kimmeta |>
  mutate(DC = ifelse(hpca_main == 'DC', 'DC', 'na')) |>
  dplyr::count(orig.ident, DC) |>
  group_by(DC) |>
  calc_frac_conf_on_grouped_count() |>
  filter(DC == 'DC')

kim_dc |> write_csv('mission/FPP/xiangya/kim_dc_count.csv')

sobj_epi <- sobj |>
  filter(tissue == 'epidermis')

# select only epidermis ----
sobj_epi |> write_rds('mission/FPP/human_SLE_epi.rds')

# find trpv3h in kim2023 -----------
sobj_kkera <- sobj |>
  filter(dataset == 'kim2023' & hpca_main == 'Keratinocytes')

sobj_kkera <- sobj_kkera |>
  quick_process_seurat(skip_norm = T)

sobj_kkera |> DotPlot('TRPV3')

sobj_kkera |> FindAllMarkers(features = 'TRPV3',
                             only.pos = T)

sobj_kkera |> summarise(n(), .by = seurat_clusters)

# CCR6 in cDC2? --------
sobj.epi <- read_rds('mission/FPP/xiangya_sle_scRNA/data/human_SLE_epi.rds')

sle.dc <- sobj.epi |> filter(bill_fine == 'DC')
sle.dc |> FeaturePlot('CCR6')

sle.dc$orig.ident |> table()

sle.dc <- sle.dc |>
  quick_process_seurat(skip_norm = T)

sle.dc |>
  DotPlot(c('CD1C','FCER1A',
            'SIRPA','CCR6'), cluster.idents = T)

sle.dc <- sle.dc |>
  mutate(dc.subtype = ifelse(seurat_clusters %in% c(11,5,7), 'cDC2',
                             'other DC')) 

sle.dc |>
  DimPlot(group.by = 'dc.subtype')

sle.dc |>
  VlnPlot('CCR6', group.by = 'dc.subtype')

sle.dc |>
  FindMarkers(features = 'CCR6', group.by = 'dc.subtype', ident.1 = 'cDC2')
